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Abstract 

The structural properties of an economical model for a confined plasma turbulence governor 
are investigated through bifurcation and stability analyses. A close relationship is demonstrated 
between the underlying bifurcation framework of the model and typical behavior associated with 
low- to high-confinement transitions such as shear flow stabilization of turbulence and oscillatory 
collective action. In particular, the analysis evinces two types of discontinuous transition that are 
qualitatively distinct. One involves classical hysteresis, governed by viscous dissipation. The other 
is intrinsically oscillatory and non-hysteretic, and thus provides a model for the so-called dithering 
transitions that are frequently observed. This metamorphosis, or transformation, of the system 
dynamics is an important late side-effect of symmetry-breaking, which manifests as an unusual 
non-symmetric transcritical bifurcation induced by a significant shear flow drive. 
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I. 



Fusion plasmas, and possibly other quasi two-dimensional fluid systems, may undergo a 
more-or-less dramatic transition from a low to a high confinement state (the L-H transi- 
tion) as the power input is increased, with the desirable outcome that particle and energy 
confinement is greatly improved due to localized transport reduction 0. In this work we 
report on a bifurcation and stability probe of an economical model for L-H transition dy- 
namics that uncovers a mechanism by which a radical change, or metamorphosis, may occur 
in the qualitative nature of the dynamics. We apply the results of this analysis to clarify 
the relationship between the structure of the model and the physics of the process that it 
describes, and draw comparisons with characteristics of L-H transitions observed in various 
experiments. 

Since 1988 there has been much progress in developing low-dimensional (low-order or 
reduced) descriptions of L-H transition dynamics and associated oscillatory phenomena 
(see, for example, Refs |, |, |, |, |, |7|, |, |, 0, 0, [T|, [U], HD, the driving force being the 
potential power of a unified, low-dimensional model as a predictive tool for the design and 
control of confinement states. For example, a model that speaks of the shape and extent 
of hysteresis in the L-H transition would help engineers who are interested in controlling 
access to H-mode. Given the many variables and parameters that could be varied around 
a hysteretic regime, it would be cheaper — i.e., save hundreds of cpu hours and/ or many 
expensive diagnostics — to know in advance which ones actually do affect the hysteresis, and 
which do not. 

To help construe the context in which low- dimensional descriptions of plasma dynamics 
are sought, it is appropriate at this stage to make some general remarks. It makes sense to 
try to find the simplest description of an evolving system that is consistent with the time 
and space scales on which one is interested in making experimental observations of that 
system. One would like the description to incorporate the qualitative nature of the system 
structure and dynamics, so that it can be used for design and control purposes and make 
useful predictions. A truly useful description usually turns out to be a low-dimensional 
system of coupled ordinary differential equations. Such descriptions are powerful because 
they are supported by well-developed theories that give qualitative and global insight, such as 
bifurcation, stability, and symmetry theory |TR |T7 |. In principle we can map analytically the 
bifurcation structure of the entire state and parameter space of a low- dimensional dynamical 
system, but this is not yet possible for an infinite-dimensional system. 

However, the quest for a low- dimensional state space that captures the qualitative dy- 



namics of L-H transitions has been problematic. It has been shown |L5], [18[ that some of 
the models cited above do not reflect salient features of L-H transitions such as shear flow 
suppression of turbulence, or are incomplete, or show profound structural discrepancies, 
although it is intuitively reasonable to expect that manifestly different models should be 
equivalent at some deeper level if they describe the same phenomena. 

By economical, or minimal, model we mean the smallest, functionally simplest, and 
mathematically consistent model that captures qualitatively the dynamical traits that are 
typically observed over many experiments in different machines. The strength and potency of 
a minimal model is just this universality; its apparent disregard for details, numbers and unit 
dimensions is sometimes perceived — wrongly — as a weakness. In keeping with this ideology 
we introduce here a consensus dynamical model that is economical in terms of variables and 
parameters, and incorporates the smallest number of rate processes of simplest functional 
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FIG. 1: Coupled rates and feedback processes that contribute to the dynamics of L-H transitions. 
Solid arrows indicate generation rates, wavy arrows dissipation, dashed arrows feedback on rate 
coefficients; black diamond indicates negative feedback. 

form needed to reflect the universally observed dynamics. If the model is successful we 
expect additional terms to have only quantitative, not qualitative or structural, effects. We 
should also be able to identify easily its limits of validity, or where it breaks down and why. 

In section II we introduce the plasma turbulence governor as a useful schema to con- 
ceptualize and represent the major contributing rate and feedback processes, relating these 
to the corresponding dynamical system. Bifurcation and stability analyses and interpretive 
discussions, with reference to reported experiments, are given in the remaining sections. In 



section |T| we begin by determining the two highest order (most degenerate) singularities in 
the system, or organizing centers. Section [IV] describes the generic bifurcation diagram and 
discusses the hysteresis and limit cycles in the system. In section [V] we illustrate and dis- 
cuss the useful properties of the two-parameter bifurcation diagram. This discussion leads 
in to section [VI], in which we determine explicitly the transcritical metamorphosis to an 
oscillatory, non-hysteretic regime. A short summary is given in section |V11| . The Appendix 
contains a derivation of the dynamical equations. 



II. 

The schematic in Fig. [I] is a primitive of a plasma turbulence governor. (The name 
is intended to refer to archetypical mechanical exemplars of feedback controllers such as 
James Watt's 1788 steam-engine governor. In |T^] a comparable scheme was called the 
"barotropic governor", in the context of quasi two-dimensional atmospheric flows.) A power 
input q creates a pressure gradient P from which the turbulent density fluctuation intensity 
TV grows at a rate with coefficient 7. The turbulence feeds energy into the poloidal shear 
flow v via the Reynolds stress a. The shear flow is generated externally at rate (p and 
damped by the ion viscosity /i. The turbulence is damped quadratically with coefficient [3. 
Also indicated is a competitive distribution of energy from the pressure gradient, whereby 
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different fractions may partition into turbulence generation and shear flow damping. It is 
not difficult to appreciate how the various rate and competitive processes in Fig. [l| could 
balance out — or rather, tin-balance out — so as to give rise to the oscillatory and hysteretic 
dynamics that are characteristic of L-H transitions. 

The reduced dynamical system that models this scheme is based on the Sugama-Horton 
model pj, which itself was derived from approximate resistive MHD vorticity and pressure 



convection equations 20, 21 



ef^-yPN (1) 

dN 
~dT 



jPN - av 2 N - f3N 2 (2) 



2— =avN - //(P, N)v + ip (3) 

//(P, N) = bP m + aP r N. (4) 
In terms of the shear flow kinetic energy F = v 2 Eqs |^ and |3] may be written as 

_ =1 PN - aFN - (3N 2 (g) 

=aFN - /i(P, N)F + ipF 1 / 2 . (g) 

The derivation of this system is given in the Appendix. The most important modification 
to the original Sugama-Horton model is the symmetry-breaking term <p in Eq. |3|. It will be 
seen that this term, which may be interpreted as an external shear flow driving rate, has 
dramatic effects on the bifurcation structure of the system. 

The first and second terms in the bipartite viscosity function, Eq. |4], model the neoclas- 
sical and anomalous viscosity damping respectively. In a plasma of low collisionality the 
exponent m is negative so a high pressure gradient has the effect of blocking the neoclassical 
contribution. (Refer to Fig. [I].) Under these circumstances energy can accumulate in the 
shear flow then feed back into turbulence decorrelation. On the other hand, a high pres- 
sure gradient and high turbulence levels both enhance the anomalous viscosity damping, 
because the exponent r is positive. The net effect will depend on the relativity of the three 
competitive rates involved in the distribution of energy from the pressure gradient. 



III. 

Generally in bifurcation analysis we are interested in the multiplicity, stability, sin- 
gularity, and parameter dependence of zero solutions of a bifurcation equation g = 
G(x, Ai, A2, • • • , A n ), where x is a state variable and the Aj are parameters, that is deriv- 
able (in principle if not always in practice) from the equilibria of a dynamical system. 

In Eqs [l]-44| we may select x = P and the principal bifurcation parameter Ai = q and set 
the right hand sides to zero to obtain the bifurcation equation, 

(apr q - qa + 6P 1+ - 7 ) (q(3 - PV) + V{P y*~?® 1/2 (5) 



2PW V * * ' ; 2(P« 7 ) 1/2 
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(where Eq. ^ has been used). Singular points occur where g = g P = 0. (Subscripts on g 
denote partial derivatives with respect to the subscripted variable.) On the v = branch 
they are given by 

{P,q,<p) = (P i ,P?'f//3,0), 

with the Pi given by the real, positive roots of f3b + P 1 ~ m (aP r — a) r y = 0. At these points 
g q = and gpp = —8 (aP[ ( — 1 — m + r) + (1 + m) a) 7 2 / (a/3) . Thus for some values of the 
exponents m and r one or more of the singularities may comply with the pitchfork conditions 

9 = 9p = 9 Q = = g PP = 0, g PPP ^ 0, g Pq ^ 0. (6) 

Obviously (since g PP must equal 0), compliance with these conditions also implies the exis- 
tence of hysteresis. 

To specify the dependence of the viscosity damping on the pressure gradient in Eq. |4] we 
set m = —3/2 and r = 1, as in [J7|. This value of m applies for the temperature dependence 
of the ion viscosity in a low collisional regime |22|| . The value of r = 1 is the simplest that is 
consistent with the suggested dependence of the anomalous viscosity on the ion temperature 
23|. 



in 



With this specification the conditions in Eq. |5] applied to Eq. |5] find the unique pitch- 
fork P* as 

/ a 2 7 2 2a 3/ y\/a/a \ 

(v, q, b, ip) = 0, — ^j, ^ ' , . P* 

^ ^ \ 9a 2 /3' 27v/3a 2 /3 J 1 ' 

At P* the two non- degeneracy conditions in Eq. || evaluate as g Pq = 8a/a, g ppp = 
— 18a7 2 /(a/?). A pitchfork is described as a codimension 2 singularity, because its uni- 
versal unfolding requires 2 parameters additional to the principal bifurcation parameter. 
Note that the second unfolding parameter, chosen here as b, can be any of the dissipative 
parameters a, b, or (3. For reference the bifurcation diagram in Fig. [| has been computed 
and plotted for the critical set (P*). The singular point on the v — branch at high q 
complies with the conditions 

9 = 9p = 9 q = 0, 9pp ^ 0, det d 2 g < 0, (7) 
where d 2 q is the Hessian matrix of second partial derivatives f ^ PP ® qP | . A singular point 

\9qP 9 qq J 

that satisfies these conditions is usually termed a transcritical bifurcation, or sometimes a 
"simple bifurcation" . 

For non-critical values of b (i.e., b ^ bn>*} = 18.58 . . . ), P* collapses to a second trans- 
critical bifurcation on the v — branch. These two transcriticals coalesce and annihilate 
each other at a second codimension 2 singularity D* on the v — branch, defined by the 
conditions 

9 = 9p = 9g = det d 2 g = 0, g PP ^ 0, g Pq ^ 0, (8) 
and found using Eq. ^| as 

v , qM= l0 lW -, _ ,0|, (D*) 
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FIG. 2: Bifurcation diagram for the critical set (P*), <p = 0, b = 18.58, a = 2.4, (3 
a = 0.3, e = 1.5. 
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FIG. 3: Bifurcation diagram with b = 1, <p = 0.02, other parameters as for Fig. |2[ 

with g PP = — 327 2 / (7/5), gp,j = 16a/ (5a). The bifurcation diagram showing this point (at 
(v, q, b, ip) = (0, 0.61 . . . , 53.52 ... , 0) for values of the other parameters as in Fig. ^) would 
be extremely dull and flat — it consists only of the line v — 0. Such a highly dissipative 
system has no interesting behaviour at all. 



IV. 



A bifurcation diagram for non-critical values of the unfolding parameters b and if is 
shown in Fig. (In the bifurcation diagrams stable solution branches are indicated by solid 
lines, unstable branches by dashed lines, and the dotted lines trace out the maximum and 
minimum amplitude of limit cycle branches.) The symmetry evident in Fig. [| is broken by 
selection of a small positive value of (p, which determines a preferred direction of the poloidal 
shear flow. 

Branches of stable limit cycles emanate from Hopf bifurcations on the +v and —v H-mode 
branches. They reflect reports from experiments that a transition to a quiescent H-mode can 
be achieved followed by the onset of oscillatory behaviour, or edge-localised modes (ELMs), 
as the power input continues to be increased 123, p6|. 27]. The original Sugama-Horton 



model was found to exhibit a chaotic time series for a particular set of parameter values in 
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FIG. 4: A time series for q = 4 on the +v branch. Other parameters are the same as in Fig. 



this regime |pq| . In our model we have found that this branch of limit cycles can undergo 
several successive period doublings followed by period halvings back to a period-one limit 
cycle. 

The limit cycles are also extinguished at Hopf bifurcations. In [HJ it was shown that 
oscillatory behaviour is regulated by the contribution of the pressure gradient evolution. At 
moderately high power input the pressure gradient and turbulence are high, neoclassical 
viscous damping is inhibited, and large amplitude oscillations would be expected as energy 
alternately accumulates in the shear flow and is exchanged with the turbulence. The relative 
phases of v, N, and P are shown in the time series of Fig. |j. However, this is balanced by 
the enhancement of anomalous viscosity damping by the larger amounts of turbulence and 
pressure gradient energy at higher power inputs. (Refer to the governor schematized in Fig. 
[I].) As this anomalous viscosity effect begins to take over the amplitude of the limit cycles 
decreases rapidly until they are extinguished at the Hopf bifurcations at higher q. 

Although definitive experiments have not yet been performed that measure the growth 
and extent of the H-mode oscillations over the power input, it is physically reasonable that 
they would be limited by some damping factor. The passage through an oscillatory regime 
with increasing power is a characteristic of type III ELMs 



29, Bflj. However, the quantitative 



features of type III ELMs, such as the frequency spectrum, are not reproduced by this simple 
model. 

On the v < branch the limit cycles are smaller in amplitude and occur over a smaller 
range of the power input. At q « 2.5, for example, the +v H-mode is oscillatory but the 
—v H-mode is quiescent. Again, to our knowledge the appropriate experiments have not 
yet been carried out, but this is reminiscent of the prescription given in ref. ||31||: "The 



key factors in creating the quiescent H-mode operation are neutral beam injection in the 
direction opposite to the plasma current (counterinjection) plus cryopumping to reduce the 
density." 

Reports of reversals in the direction of main or impurity ion poloidal shear flow E3, 15B 
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FIG. 5: Bifurcation diagram with N as chosen state variable, showing the curve that corresponds 
to the positive v branch in Fig. ^. 



can also be rationalized on the basis of Fig. |3]. In a system that is evolved initially onto the 
v < branch, the poloidal shear flow must reverse if a perturbation decreases the power input 
slightly below that at the lower limit point. Shear flow reversal may also occur anywhere 
along the v < branch, if the system is given a sufficiently strong transient kick. 

Note that in Fig. ||| the shear flow v reaches a broad maximum with increasing power 
input, then decreases to the the pre-transition level given by the shear flow source. This 
would be reasonable behaviour on physical grounds — one would not expect the shear flow to 
increase indefinitely with power input, because the turbulent viscosity damping (the second 
term in Eq. |] with r = 1) begins to take over as the power input increases the pressure 
gradient. 

Clearly there is scope for tuning other parameters in the model so as to obtain a complete 
picture of the steady states and limit cycles over parameter space, and more quantitative 
agreements with experiments. One may wish, for example, to maximize the range of q over 
which turbulence stabilization occurs, or minimize the range of q over which limit cycles 
occur, or both. 

Figure [5] also shows the hysteresis that is predicted by compliance with the conditions 



in Eq. H. Transitions with hysteresis have been observed in several machines: DIII-D p5 



Asdex Upgrade |34], pSfl , JET, and in simulations of ITER p6fl , and Alcator C-Mod [|36 
Hysteresis is typically modified by dissipation, characterised in this model by the parameters 
/3, b, and a. However, hysteresis does not seem to be a necessary or universal feature of 
discontinuous transitions. 

One of the typical features of L-H transitions that a minimal model should reflect is 
suppression of the turbulence by the shear flow. Figure ^] shows the bifurcation diagram 
with the mean square turbulence level N as the state variable, where for clarity only the 
curve that matches the positive v branch is given. The turbulence is clearly suppressed over 
the hysteretic region, then begins to grow again as the higher pressure gradient from higher 
power input creates more turbulence. 



V. 



The width and extent of hysteresis for selected values of b can be judged from the two- 
parameter bifurcation diagram for the +v branch in Fig. in which computed curves of 
the singular points in Fig. ^| are shown. The solid lines mark the loci of limit points (which 
are also sometimes called fold or saddle-node bifurcations) as b is varied. The dot-dash line 
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FIG. 6: Two-parameter curves of the singular points in Fig. 
points, dot-dash lines are the Hopf bifurcation loci. 
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Solid lines are the loci of the limit 



is the locus of Hopf bifurcations as b is varied. If one can imagine taking slices across this 
diagram at various important values of b, the bifurcation story of the system can be told 
compactly, by inferring a reconstruction of the single-parameter (q, v) bifurcation diagram 
corresponding to each selected value of b. 

A slice taken above the critical value of b at the cusp would yield a (q, v) bifurcation dia- 
gram that shows no multiplicity of states. Thus, in a highly dissipative system any transition 
is expected to be smooth and gradual rather than discontinuous, and a number of experi- 
ments suggest this conjecture. In ASDEX Upgrade the power hysteresis disappears at higher 
density (which implies more collisional damping) where gradual rather than discontinuous 
confinement improvement occurs [05 |. A regime in which density fluctuation amplitudes 
are reduced continuously was also observed in |37 



In fl3ql a discontinuous bifurcation of 



the electric field in a stellarator was reported for conditions of low neutral density, where 
the charge-exchange damping rate is low. The change in the electric field became gradual 
for conditions of high neutral density, because the charge-exchange damping rate increases. 
(The electric field is related to the poloidal shear flow and the pressure gradient through the 
radial force balance fl39f .) 

Oscillatory behaviour is also expected to be damped out at high dissipation rates. The 
maximum in the locus of Hopf bifurcations in Fig. ^| occurs at the value of b where the 
two Hopf bifurcations on the +v branch in Fig. [5] annihilate each other (or conversely, are 
created). Above this value of b the +v branch is stable with no associated limit cycles. 

As slices are taken at lower b the hysteresis and the range of oscillatory behavior evidently 
become broader. At low dissipation rates the feedback is strong and nonlinear behavior is 
expected to be more pronounced. 

The crossing of the Hopf and limit point loci in Fig. ^| is non-local, i.e., the value of P 
(and of v and N) at the crossing on the Hopf curve is different from that on the limit point 
curve. Within the overlapping region a direct transition to an oscillatory H-mode may occur. 



VI. 

The minimum in the limit point curve of Fig. ^| implies the existence of another tran- 
scritical bifurcation (defined by Eq. |^) in the system, that occurs at non-zero v. This 
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FIG. 7: A series of bifurcation diagram snapshots taken at increasing values of ip illustrates the 
exchange at fT m and its aftermath. Here e = 1.0 and other parameters are the same as in Fig. 0. 



non-symmetric transcritical bifurcation may have important issues concerning access to and 
control of confinement states. Consider the series of stills in Fig. |7|, which snapshot +t> 
bifurcation diagram sections for increasing values of <p. 

As if is incremented a separate branch of solutions to Eq. ^, which is trapped at (q, v) = 
(0, oo) for ip = 0, begins to intrude more prominently into the physical region (if = 0.05). 
Although it is unstable at first, and therefore physically irrelevant for very small values of 
if, it does not remain so. The singular occurrence of a zero real eigenvalue and a pair of 
complex conjugate eigenvalues with zero real part signals the appearance of a degenerate 
Hopf bifurcation. (Eigenvalues were computed numerically.) 

Further increments in if separate the limit point and the Hopf bifurcation, between which 
the solutions are stable (if = 0.08). The branch of limit cycles that emanates from the Hopf 
bifurcation undergoes one or more period doubling bifurcations before ending, presumably 
at a homoclinic (infinite period) terminus. This branch of limit cycles is quite short, and 
thus not very well resolved in Fig. |7| for the lower values of (p. Note also that the branch of 
limit cycles emanating from the hysteretic solution branch has also appeared by (p = 0.08. 

At a metamorphic value of ip that we designate ifTm the arms of the two separate 
steady-state branches are exchanged at a transcritical bifurcation. We know this point 
is present because the denning conditions Eq. [7] are satisfied with <p ^ 0. (In numbers 
(v,q,<f) Tm = (1.8247.., 0.1468.., 0.08059..), with det d 2 g = -0.004687.., g PP = -0.001250.., 
for values of the other parameters as given in Fig [| Note that the value of e is irrelevant 
for calculating steady-state bifurcations such as the pitchfork and transcritical but not for 
Hopf bifurcations.) 

After the exchange, e.g. at ip = 0.085 and ip = 0.1, we see the unusual occurrence of three 
Hopf bifurcations on the same branch, although this situation could be inferred from the 
shallow but distinct minimum and the maximum in the locus of Hopf bifurcations in Fig. |6|. 

The last frame of Fig. |7], for if = 0.11, is taken after the "new" and the lower-g "old" 
Hopf bifurcations have collided and annihilated each other at a singular point associated 
with two zero eigenvalues. This is what the minimum in the curve of Hopf bifurcations 
means. The branch of limit cycles emanating from the upper-g "old" Hopf bifurcation now 
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continues to the (presumed) homoclinic terminus. There are a couple of period-doublings on 
it (not shown). We also see that the limit cycles are extinguished and a quiescent H-mode 
is achieved at the single remaining Hopf bifurcation. 

Turning our attention to the stable part of the lower branch in the last frame of Fig. [7], 
we see that as q is tuned past the lower limit point the system must jump to another 
stable attractor. This transition is very different from the intrinsically hysteretic transition 
depicted in Fig |3]. Here the stable attractor on the upper branch is a limit cycle rather than 
a fixed point. Furthermore, this transition is not hysteretic. In fact, hysteresis is (locally) 
forbidden by the condition gpp ^ of Eq. [7|. Therefore it is not modulated by dissipation in 
the same way as the transition in Fig.§, although the feedback itself is still due to nonlinear 
dissipation rates. 

As the value of <p is increased even further, bifurcation diagrams that one could plot 
gradually become less meaningful. This is because ip = constant is a first approximation, 
valid for small (p, to to a nonlinear function <p(C), where £ may include dynamical variables 
and parameters. 

This type of transition could serve as a model for the dithering or L-H-L transitions, 
followed by a quiescent H-mode, that have been reported in many machines. Although there 
may be other mechanisms for dithering transitions — another possible scenario is given at the 
end of section [V| and indicated in Fig. || where an oscillatory transition may occur in a very 
poorly dissipative system — we have at least a preliminary semiological and classification 
guide: if your transition is oscillatory and non-hysteretic then perhaps you should look for 
a strong shear flow source, if it is strongly hysteretic perhaps you should look at dissipation 
mechanisms. Some experimental evidence supports the idea that dithering transitions result 
from a strong shear flow source. In [ i0|] , an analysis of time series data around the L-H 



transition in COMPASS-D suggested that a homoclinic orbit is involved in the change of 
stability at the transition. In stellarator W7-AS typically the quiescent ELM-free H-mode 
is obtained after a phase characterized by quasi-periodic ELMs [il, i2|. In HI stellarator 
a transition to fluctuating H-mode occurs at lower gas filling pressures and lower magnetic 
fields than the transition to quiescent H-mode |43|| . 

In terms of the governor in Fig. [I] a shear flow that is generated internally and driven ex- 
ternally at comparable rates is likely to give rise to interesting non-linear dynamics, because 
more kinetic energy in the shear flow leads to more turbulence suppression through decor- 
relation, but also a larger damping effect, which then alters the competitive distribution of 
energy from the pressure gradient. 



VII. 

In summary, this reduced dynamical model, comprised simply of energy input, exchange, 
and loss rates, reflects generic characteristics of confined plasma bulk dynamics that have 
not been reflected in previous models. The bifurcation and stability analysis also reveals 
two qualitatively different transitions. The hysteretic transition is controlled by the damping 
rate coefficients. The non-hysteretic transition occurs when there is a relatively strong shear 
flow drive. 

Symmetry-breaking in this system has two major effects. Firstly, a non-zero shear flow 
drive is physically inevitable, even in the best-controlled experiments, and it determines a 
preferred direction for the shear flow. Secondly, it interacts with the internal generation and 
loss dynamics to cause the metamorphosis shown in Fig. [7|. 

More generally, the information obtained from this analysis strengthens the thesis devel- 
oped in [T7[] : that remarkably low-dimensional models can capture and help explain essential 
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aspects of turbulent flows that elude understanding from numerical simulations that include 
resolved spatial scales, and that physical deductions can be made from observations of bi- 
furcations. 
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FIG. 8: Simple slab geometry is assumed. The plasma edge region is —5 < x < 5, with x = 5 at 
the plasma surface. Vpo < is the y, z-averaged pressure gradient. 



Appendix 

Reduced MHD fluid equations in tokamak and stellarator geometries were originally de- 
rived by Strauss |20], 21] . In the electrostatic approximation, a damped MHD fluid may be 
described by the following momentum and pressure convection equations: 



p^- = -Vp + J x B + fiVjV + Vt'px - pu(v -V (x) y) (9) 

|=xVip, (10) 

where d/dt = d/dt + v • V, together with the incompressibility condition V ■ v = and 
the resistive Ohm's law E + v x B = 77J. The symbols and notation are explained in 
table |. The curl of Eq. ^ yields a vorticity equation, which is sometimes preferred in 
two-dimensional fluid dynamics, but we have used the momentum form because it is more 
transparent physically and has a simpler correspondence to the kinetic energy. An infinite 



slab configuration is used for simplicity and generality, as was also assumed in |^3| for a drift- 
kinetic treatment of plasma relaxation. It is sketched in Fig. |8], where the region —5 < x < 5 
can be taken to represent a region at the edge or within a confined plasma where a transport 
barrier evolves. 

The last term on the right hand side of Eq. |S] removes the nonlinear shear-flow rever- 
sal symmetry of the system under v x (x,y,t) — > v x (x, — y, t), v y (x,y,t) — > —v y (x,—y,t), 
p(x,y,t) — > p(x,—y,t), V(x) — > V(x). Similar equations, without the symmetry-breaking 
term in the momentum balance, have been used by several authors [0, |23|, |47] 



cLS cL 



basis for studying resistive turbulence-flow interactions. The symmetry-breaking term was 
introduced in [p3[| , but only a posteriori as an adjunct in an equation for the background 
poloidal flow. Here we introduce it at the outset. It models the friction force acting between 
the single-fluid plasma velocity v and an assumed external poloidal flow Vy. Although Vy 
may be described for convenience as an external velocity, the term represents any asymmet- 
ric shear-inducing mechanism, such as friction with neutrals, non-ambipolar ion orbit losses, 
or neoclassical effects not included in the slab model. 

The symmetry operation on t>o(x) = (v y )(x) and V(x) is sketched in Fig. |[ We are 
working in the frame in which there is no electrostatic potential difference across the slab. 
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FIG. 9: Without the friction force the system is invariant under the transformation vo(x,t) — > 
—vq(x, t) (solid line), V(x) — > V(x) (dashed line). When the friction coefficient v ^ the symmetry 
is broken. 



That is, it is assumed that we have made a Galilean transformation to the frame in which 
the spatial average of Vq across the slab is zero. For simplicity we also assume that the 
spatial average of V is zero. 

Equations and ETUI are not intended to express a detailed fluid description of a plasma, 



but are intended instead to represent a qualitatively authentic, semi-empirical model for the 
essential generation and loss processes that give rise to the turbulence-shear flow interactions 
that we have schematized in Fig. [I] as the plasma turbulence governor. The dynamical 
system Eqs can be derived from Eqs || and [LC], following the spatial averaging procedure 
implicit in J7|. 

First of all, the dynamics of the mean flow vq = (v y ) are extracted from the first moment 
(v y y) of ( Eq. |) as 



d t v - fid^vo + d x (v x v y ) = -u(v - V), 



(11) 



the energy moment of which gives the spatially averaged evolution of shear flow kinetic 
energy F: 



d 
dt 



This may be written as 



dx 



dF 
~dt 



dx 



+ 



_ dv , 



dx 



dx (v x v y ) 



-s 



dvo 
dx 



(12) 



+ - / dx vVvq. 

.1-8 



-e F + E F + E ( 



(13) 
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where the definitions of ep, Ep, and £L correspond respectively to each term on the right 
hand side of Eq. and F = ± J S _ S f vg. 

Next, the second moment of Eq. |9] gives the total rate of evolution of F and turbulent 
kinetic energy N: 



d 

dt. 




(14) 



+ T / fife ^^Wo, 

<5 ./-a 



which may be expressed succinctly as 



d 



— [F + N] = E N - e N - e F + E v 



(15) 



where E N and are defined by the first two terms on the right hand side of Eq. [TJ] and 

iV — 8 J-S 2 17 ' 

Finally, the evolution of potential energy in the pressure gradient is defined from the 
ic-moment of Eq. [H], assuming the cross-field thermal transport can be neglected: 



d 
dt 



1 



dx (—x) Q 



,Po 



-s 



5 dx <&La, 

-5 P 



or 



dP 
~dt 



Ep-E 



N ■ 



(16) 



(17) 



with the input rate Ep defined as the first term on the right hand side and P = 
tftgdx {-x)Q'p /p. 

The spatially averaged dynamical system thus consists of Eqs [15], 0, and [I7|. For closure 
we follow [J/J, using the approximations po(x) ~ p (x = 5) + xdp /dx and v (x) ~ v (x = 
5) + xdvo/dx for the background pressure and flow profiles and re-defining P and F as the 
gradient terms alone. Approximations or expressions based on empirical arguments were 
given in for the rates in Eqs 0, |lf|, and 0. The rates given in Eqs [l|-|3| are economized 
versions of those expressions, in the sense that simpler power laws were chosen if this did not 
result in any qualitative changes to the singularity and stability structure of the system. The 
rationale is that for most of the rates we shall only learn from experiments whether different 
powers apply, meanwhile simple power laws give more transparent algebra. We approximate 
the energy transfer rate from the pressure gradient simply as E^ ~ ( / ~f/e)PN, and the 
energy transfer rate between the turbulence and the shear flow, due to the Reynolds stress, 
as Ep ~ aFN. The power input through the boundary is defined as E p = q/e. The two- 
timing coefficient e is related to the thermal capacitance, and regulates the contribution of 
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the pressure gradient to the dynamics. For the dissipative terms we take the turbulent energy 
dissipation rate as ~ f3N 2 and the shear flow energy damping rate as ep — fi(P,N)F, 
assuming the viscous damping to be dominant in ep. The external shear flow driving rate 
is then E v ~ pF 1 / 2 , with ip ~ bvV . To obtain the evolution of the shear flow in terms of a 
velocity gradient variable we re- define v = ±F 1 / 2 . Eqs |-@ ensue. 

v = -g-z x V</> = vq + v E x B flow velocity 



v o = ( v ) average background component 

v fluctuating or turbulent component 

P = Po + P plasma pressure 

Po = (p) average background component 

p fluctuating or turbulent component 

p average mass density of ions, assumed constant 

fi ion viscosity coefficient 

Bo magnetic field along the z axis 

n resistivity 

v frictional damping coefficient 

SI' = dfl/dx > average field line curvature, assumed constant 

vi d 2 + d 2 

V|| d Z + f-dy 

X cross-field thermal transport coefficient 

V external flow 

(• • • ) average on (y, z) plane 



TABLE I: Glossary of symbols, terms, and notation. 



